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A solution to an equilibrium of irrotational binary polytropic stars in Newtonian gravity is expanded 
in a power oi e = aa/R, where R and ao are the separation of the binary system and the radius of 
each star for R — oo. For the polytropic index n = 1, the solutions are given almost analytically 
up to order e^. We have found that in general an equilibrium solution should have the velocity 
component along the orbital axis and that the central density should decrease when R decreases. 
Our almost analytic solutions can be used to check the validity of numerical solutions. 
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Coalescing binary neutron stars (BNSs) are consid- 
ered to be one of the most promising sources of gravita- 
tional waves for laser interferometers such as TAMA300, 
GEO600, VIRGO and LIGO We can determine the 
mass and the spin of neutron stars from the gravitational 
wave signals in the inspiraling phase. We may also ex- 
tract the informations on the equation of state of a neu- 
tron star from the signals in the pre- merging phase 
For this purpose it is important to complete theoretical 
templates of gravitational waves in the pre-merging phase 
as well as in the inspiraling phase. 

Recently Bonazzola, Gourgoulhon and Marck have nu- 
merically calculated quasi-equilibrium configurations of 
irrotational BNSs in general relativity Q^. Although 
their results seem to be reasonable, we do not have any 
analytic solutions in general relativity in order to check 
the validity of their results. 

On the other hand, Uryii and Eriguchi have numer- 
ically constructed stationary structures of irrotational 
BNSs in Newtonian gravity . We can use semi-analytic 
solutions produced by Lai, Rasio and Shapiro (hereafter 
LRS) for checking the validity of the results. 

However, in numerical solutions of Uryti and Eriguchi, 
the velocity component along the orbital axis exists while 
in those of LRS such a component is assumed to be zero 
from the beginning. When we extend the analytic solu- 
tions to the general relativistic ones (D , we should include 
this velocity component. This is because in the numer- 
ical calculation, there is a possibility to obtain another 
solution although the binding energy of a BNS is almost 
the same value, and to lead a different conclusion 

In order to include the velocity component along the 



*The irrotational state is considered to be a realistic state 
for binary neutron stars before merger Q,^. 
^Present address: Departement d'Astrophysique Relativiste 
et de Cosmologie, UPR 8629 du C.N.R.S., Obserbatoire de 
Paris, F-92195 Meudon Cedex, France 



orbital axis, we solve the equation of continuity with the 
other basic equations. The method we use in this Letter 
is that we seek a solution to an equilibrium of irrotational 
binary polytropic stars in Newtonian gravity by expand- 
ing all physical quantities in a power of e = ao/i?, where 
R and ao are the separation of the binary system and the 
radius of each star for R — oo. We extend the method 
developed by Ghandrasekhar more than 65 years ago for 
corotating fluids to the one for irrotational fluids. 

Although a binary system consists of two stars, we pay 
particular attention to one of two stars. We call it star 1 
whose mass is Mi and the companion one star 2 whose 
mass is M2. In this Letter, we adopt two corotating co- 
ordinate systems. First one is X whose origin is located 
at the center of mass of the binary system. For calcula- 
tional convenience, we choose the orbital axis as X3, and 
we take the direction of Xi from the center of mass of 
star 2 to that of star 1. The second coordinate system is 
the spherical one r = (r, 9, ip) whose origin is located at 
the center of mass of star 1. We use units of G = 1. 

Since we treat irrotational fluids in Newtonian gravity, 
the basic equations are the equation of state, the Eu- 
Icr equation, the equation of continuity and the Poisson 
equation: 
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where P, p, n, U and are the pressure, the density, 
the polytropic index, the gravitational potential and the 
orbital angular velocity, respectively. The gravitational 
potential U is separated into two parts, i.e., the contri- 
bution from star 1 to itself C/^~*^ and that from star 2 to 
star 1 C/^^^. U'^'^^ is written as 
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Here, Pi and denote the Legendre function and the 
reduced quadrupole moment, and the superscript ' means 
the term concerned with star 2. ii' is a constant related to 
entropy and v represents the velocity field in the inertial 
frame. In the irrotational fluid case, we can express v as 
a gradient of a scalar function $, i.e., v = V$. 

Following the Lane-Emden equation we first express 
the density as p = Pc©" with pc being the central density. 
By using a = [if (l+n)/9c^"~"^/(47r)]-^/^, we also introduce 
^ = r / a. We expand in a power series of a parameter e 
as 9 = X^So f*©i- Since the shape of star 1 is spherical 
when R is large, the lowest order term of O is the solu- 
tion of the Lane-Emden equation. Then we expand Qi 
by spherical harmonics as 9^ = ,^ ^^VimlO^i™!^' 
The radius of a spherical star ao is given by oq = a^i as 
usual. Note here that = tt for n = \ case. 

Now we consider the orbital motion of star 1. In the 
spherical coordinate system, it becomes 

n^r = nR{nxi)^^^ + na{nxi)^^^, (6) 

where 

X = — (sin 6* sin (/?, cos0sin(y9, cos(^), (7) 

(nx^)^,^^ = (0, 0, Csin^), (8) 

and p = Mi/M2- The first term on the right-hand side 
of Eq. comes from the orbital motion of the center of 
mass of star 1 and the second term comes from the fluid 
motion around the center of mass of star 1. 

Next, we rewrite the equation of continuity (0) as 
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The condition for $ at the stellar surface is (V'l' — Jl x r) ■ 
(V9)|^^^^= 0, since 9 = at the surface. We expand 

<& also as $ = X^i^o^*^'- gradient of the lowest 

order term of $ should agree with the orbital motion of 
the center of mass of star 1, i.e., ^IR{^1 x ^)orb because 
when R is large, the shape of star 1 is spherical and star 
1 has only the orbital motion of the center of mass in 
the inertial frame with no intrinsic spin. This leads us to 
normalize $ as <& = e<I>/(riaao). We again expand 4>i by 
spherical harmonics as l>j = Y.i,m '•*V;m(0^/'"(^> f/')- 

The orbital angular velocity is derived from the first 
tensor virial relation defined by pT| 
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where xi = rsin0cos<y9. If we substitute Eq. (||), 9o 
and $0 into Eq. (|lO|), we obtain the orbital angular 
velocity in the lowest order as JIq = Mtot/R^ where 
^tot — Ml + A/2. Note here that we also expand 

Finally, we express the gravitational potential by 
rewriting Eq. (^ as 
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U = K{l + n)p- 

where Uq is constant. Substituting Eq. ( [ll| ) into the 
Poisson equation (|4[), we obtain the equation to deter- 
mine the equilibrium figure as 

^ * ^(V$)2 - 2(V*) • (n X r) 
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Now a solution can be obtained iteratively. Firstly, 9^ 
is determined by demanding that the gravitational po- 
tential and its normal derivative are continuous at the 
stellar surface |l^, that is, Uint\^=s = Uext\^=s and 
dU.int/d£,\^=E = dUext/d£,\^=E, where S(6', (p) expresses 
the surface (9(^ = 'E.{9,(p)) — 0). Substituting 9^ and 
Eq. (|) into Eqs. (|) and (|l^, we obtain and fJ^. Af- 
ter that we substitute these equations into Eq. (|l2|) and 
derive 9i+i. We continued this procedure up to order 
in this Letter. 

Since we would like to present almost analytic results, 
we only calculate n = 1 case in this Letter. For other 
polytropic indices, we must solve differential equations 
numerically. The results of these cases will be given in 
the subsequent paper p^ . The density profile up to 
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Here ji and P™ denote the spherical Bessel function and 
the associated Legendre function. The function ^^'>'ip22 is 
obtained by solving a differential equation: 
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where ^''•'02(0 '^i^l be defined later. We point out that 
the terms Qi and 82 disappear in the equation for Q 
since there are no tidal terms to produce 61 and 62 in 
the gravitational potential. 

We can write the velocity potential as $ = $o + e''^4 + 
e'^^s + e^$6, where 
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The functions '^^2; '^'(^a and '^^^(j)A are determined by 
solving following differential equations: 
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where 0q denotes dQ^/d^. Note that the terms $1, $2 
and $3 disappear in the equation for $. In Table I, 
we show the results of the velocity potentials and their 
derivatives for an identical star binary [p — 1). From the 
expressions of $4, <&5 and $6, it is clear that the velocity 
has non-zero components along the orbital axis. Because 
since (^'02 is not proportional to , there remains the ve- 
locity component d^i/dx^ where is = ^cos0[]. We can 
also show the existence of the velocity components along 
the orbital axis for $5 and $6- This analytic result does 
not agree with the semi-analytic solution given by LRS 
but qualitatively agrees with the numerical solutions 
given by Uryii and Eriguchi ||^. 

The orbital angular velocity is calculated as 



where we assume that a black hole is a point source. 

Now it is ready to calculate relevant physical quanti- 
ties. The mass of star 1 is calculated as Mi = d^xp = 
47r/9ca^Ci[l + 45e^/(2p2^i)]. Since we consider the se- 
quence for given baryon mass, the mass of star 1 should 
be the same one, which yields 

Pc-Pco(l-2^e^), (28) 

where is the central density for R = 00. We note that 
the central density of star 2 should be written by chang- 
ing p into 1/p in Eq. (|2^). The equation ( |2^ ) means that 
the central density of star 1 up to order decreases from 
that of the spherical star as the separation decreases. Al- 
though we did not assume ellipsoidal figures like in Ref. 
0, the result is essentially the same as that derived by 
Lai ijl^, since the calculation of the central density up 
to 0(e®) includes only quadratic terms. In both calcu- 
lations the central density of irrotational binary system 
decreases as . This dependence is different from that 
of corotating ones, i.e., R^^ ||lOt| . 
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FIG. 1. The total energy as functions of the orbital separa- 
tion for an identical star binary (Afi = M2 = M). Solid line 
denotes our result. Filled circles and open triangles are the 
results of LRS (1994) and Uryu and Eriguchi (1998). 
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where = -i-n/e^. The effect of the quadrupole mo- 
ments in ft'^ is O(e^). One can see from Eq. (||) that the 
quadrupole term is 0{e^) higher than the monopole term 
since J^n is O(e^). In the case of a different mass binary, 
i.e., a neutron star-black hole binary, we can express the 
orbital angular velocity as 
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^In the case of n = 0, ^•^^4>2 is proportional to 



Next, we calculate the total energy of the binary 
system, which is written as i? = Iltot + {Wseif)tot + 
{Wint)tot+Ttot, where Iliot, {Wseif)tot, {Wint)tot and Ttot 
denote the total internal energy, the total self-gravity en- 
ergy, the total interaction energy and the total kinetic 
energy, respectively. They are given as 
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where the reduced quadrupole moment of star 1 is calcu- 
lated as = (2Afiag/3p)(15/Ci - 1). The reduced 
quadrupole moment of star 2 is expressed by chang- 
ing Ml into M2 and p into 1/p. We note that the 
above energies satisfy the virial equation for n = \; i.e. 
SHtot + {Wseif)tot + {W^nt)tot + mot = 0. Accordingly, 
the total energy is given as 




(33) 

In the case of a different mass binary, we obtain the total 
energy as 
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FIG. 2. The total energy as functions of the orbital sepa- 
ration for a different mass binary (p = 0.1). Solid line and 
open triangles are the same notations in Fig. 1. Filled circles 
are the results of LRS (1993). 

Our solution is correct if e ^ 1 so that it can be used to 
check the validity of numerical solutions. For any numer- 
ical codes, one can ask to solve an equilibrium for large 
R. One can compare numerically derived density and 
velocity distribution with our almost analytic solutions. 
However, for small i?, our expansion up to may not 
be enough. We can include the effects of the quadrupole 
terms of stars in the physical values such as the total en- 
ergy at order e^. Since the next higher order terms are 
the octupole ones at order and their coefficients are 
order unity, we think that this expansion converges for 
the effect of the deformation. However, at order e^, there 
appears the spin kinetic energy term in the total energy. 
There is a possibility to change the behavior of the to- 
tal energy for small R. Therefore, in order to apply our 
solution in this case, further higher order calculations in 
our scheme as well as the check of numerical codes us- 
ing our almost analytic solutions for large R are urgently 
necessary. 



In Figs. and H, we show the total energy for bi- 
nary systems with p = \ and p = 0.1 as functions of the 
orbital separation for an identical star binary. We find 
from these figures that the results of numerical and semi- 
analytic calculations coincide with our analytic solutions 
up to rather well in the region i?/ao > 3 for p = 1 and 
R/ao > 5 for p = 0.1 although they are quite different in 
details. 
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TABLE I. Velocity potentials and their derivatives for an 
identical star binary (p = 1). 
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0.0 


0.00 0.00 


0.00 0.00 


0.00 


0.00 


0.5 


2.46(-2) 9.86(-2) 9.64(-4) 5.80(-3) 


-3.40(-5) 


-2.73(-4) 
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9.91(-2) 0.200 


7.80(-3) 2.36(-2) 


-5.52(-4) 


-2.23(-3) 


1.5 


0.226 0.307 


2.68(-2) 5.50(-2) 


-2.87(-3) 


-7.83(-3) 
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0.408 0.424 


6.55(-2) 0.103 


-9.40(-3) 


-1.96(-2) 


2.5 


0.652 0.556 


0.133 0.171 


-2.41(-2) 


-4.13(-2) 


3.0 


0.967 0.710 


0.241 0.269 


-5.33(-2) 


-7.87(-2) 




1.07 0.760 


0.282 0.303 


-6.55(-2) 


-9.34(-2) 
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